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Abstract. We describe an extension of the cosmological hydrodynamics code Enzo to include the 
self-consistent transport of ionizing radiation modeled in the flux-limited diffusion approximation. 
A novel feature of our algorithm is a coupled implicit solution of radiation transport, ionization 
kinetics, and gas photoheating, making the timestepping for this portion of the calculation resolution 
independent. The implicit system is coupled to the explicit cosmological hydrodynamics through 
operator splitting and solved with scalable multigrid methods. We summarize the numerical method, 
present a verification test on cosmological Stromgren spheres, and then apply it to the problem of 
cosmological hydrogen reionization. 
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INTRODUCTION 

It is a distinct privilege to help celebrate Dimitri Mihalas' 70th birthday by reporting 
on our latest work at this Festschrift in his honor. As a former faculty colleague of his 
at UIUC and scientific collaborator, I could tell many anecdotes. However of the many 
things I learned from Dimitri, the two that stick with me are: (1) solve the problem you 
need to solve to do the science, and (2) be rigorous about how you do it. These lessons 
motivate our present approach to self-consistent cosmological radiation hydrodynamics 
which is documented more fully in 

A current frontier in cosmological structure formation simulations is including the 
feedback of radiating sources such as galaxies and AGN on the intergalactic medium 
(IGM) in a self-consistent way. For example, the collective UV radiation from proto- 
galaxies is believed to reionize the IGMatz^7-8^3]. This process can be thought of 
as the expansion and eventual overlap of R-type ionization fronts driven by high rates 
of star formation in the protogalaxies. R-type ionization fronts couple to the gas very 
weakly. Consequently a number of studies have simulated cosmological reionization by 
post-processing density fields taken from cosmological simulations with a standalone 
radiative transfer code ; e.g. [0, 3, [s]]. However, closer to the sources and within the 
galaxies themselves where the gas density is higher, or when an intergalactic R-type 
front sweeps over a dense clump, ionization fronts may become D-type. When this hap- 
pens coupling to gas motions is strong and a self-consistent approach to modeling is 
required in which hydrodynamics, radiative transfer, and the thermal/ionization state of 
the gas are evolved in a coupled fashion [@]. 



Ivj] summarizes a variety of methods currently under development by the numerical 
cosmology community that do this. The radiative transfer methods employed include 
ray tracing, Monte Carlo, and moment methods. Here we present a method based on 
flux-limited diffusion. A novel feature of our algorithm is a coupled implicit solution of 
radiation transport, ionization kinetics, and gas photoheating, making the timestepping 
for this portion of the calculation resolution independent. This will be essential when 
adaptive mesh refinement (AMR) is employed. At present our algorithm is only im- 
plemented on uniform (non-adaptive) Cartesian grids. After describing our method, we 
verify its correctness on a cosmological Stromgren sphere test problem. We then present 
a low-resolution "first light" test of the coupled code to the problem of cosmological 
reionization. A more comprehensive description of our method and verification tests can 
be found in [HI]. 



PHYSICAL MODEL 

We consider the coupled system of partial differential equations 

<?fpfc + -Vfc- Vpi = --pfoV-v^, (1) 
a a 

dtyb + -{yb-^)yb = --yb Vp--V(^, (2) 

a a apb a 

1 1 1 

(9,e + -Vfe-Ve = — -e V • (pv^) - -v^ ■ + G- A (3) 

a a apiy a 



1, 



dtvii + -V • {mvb) = auneTij - n/rf\ i=l,...,Ns, (4) 



a 



dtE + -V-{E\h)=V-{DVE)--E + 47tri-cKE, (5) 
a a 

y^^ = ^{Pb+pd,n-{p))- (6) 

These describe conservation of mass Q), conservation of momentum conservation 
of energy ([3]), chemical rate equations ([H), flux-limited diffusion radiative transfer dS]) 
and Poisson's equation in a coordinate system that is comoving with the expanding 
universe [l8|,@,[l^(lj. The independent variables in these equations consist of the comov- 
ing baryonic density pb, the proper peculiar baryonic velocity v/,, the total gas energy per 
unit mass e, the comoving number density for each chemical species n,-, i = 1, . . .,Ns, 
the comoving radiation energy density E, and the modified gravitational potential ^. 

The cosmological flux-limited diffusion (FLD) equation ^ deserves some comment. 
In deriving this equation from the general multi-frequency version ifioll . 

dtEy + ^V -{Evyb) = ^ ■{DVEy)+v^dvEy+47rriy-cKyEy, (7) 



we have assumed a prescribed radiation frequency spectrum, XEi^), allowing the 
frequency-dependent radiation energy density to be written in the form Ev{x,t,v) = 



E{xj) XE{y)- With this assumption, the single "grey" radiation energy density is given 
by 

/"CO /-oo 

E{x,t)= Ey{x,t,v)dv = E{x,t) XE{y)dv, (8) 

and the equation ([5]) is then obtained through integration of O over frequencies rang- 
ing from the ionization threshold of hydrogen (hVQ = 13.6 eV) to infinity. In this 
paper we assume the radiation has a Tg = 10^ blackbody spectrum, i.e. Xsiv) = 

8^M^)V(-p(&)-i). 

The dependent variables in these equations are the proper pressure p, the temperature 
T, and the comoving electron number density Hg, given through the equations 

p=pi,{r-n(e-^\vi,n, (9) 

r = (r-,)i^, (10) 

J n/f// , (hydrogen only) 

rig = < , I (11) 

[T^Hii + i^Heii + j^Heiii: (hydrogcn + hehum). 

Here 7 is the ratio of specific heats, which we take to be 5/3; corresponds to the mass 
of a proton and kj^ is Boltzmann's constant. The local molecular weight /i depends on 
the density and chemical ionization state. 

In addition, the system ([B-® contains a number of coupling terms and coefficients. 
The coefficient a{t) = (1 + z)^^ denotes the cosmological expansion parameter for a 
smooth homogeneous background, where the redshift z is a function of time only; all 
spatial derivatives are therefore taken with respect to the comoving position x = r/a{t). 

The term d = a{t) is obtained by integrating the Friedmann equation for the set 
of assumed cosmological parameters. The gas heating and cooling rates G and A are 
functions of the temperature, radiation energy density and chemical ionization state. 
The temperature-dependent chemical reaction rates define the interactions between 

chemical species, and the photoionization rate Tf'^ depends on the radiation energy 



density. Formulas for all of these terms may be found in the references [|lll . I12L 113 . 

In the radiation equation ©, c is the speed of light, the total opacity K is a function 
of the chemical ionization state, and the emissivity rj is provided as either a radiation 
source term, or may depend on the density, velocity, gas energy, and chemical ionization 



state. The formulae for these dependencies may be found in the references 11121 11711. Of 
special importance in this equation is the coefficient function D, which in a flux-limited 
diffusion approximation attempts to allow the equation to span behaviors ranging from 
nearly isotropic to free- streaming radiation. To this end, we choose the coefficient to be 
of the form 

D(£)=diag(Di(£),D2(£),Z)3(£)), where A(£) = -^l^^^, (12) 

OKf + 3KTRi+Rf 



with Rj = \djE\/E,i = 1,2,3. Here Kj = K+ Ks is the total extinction coefficient, 
where k is the opacity and Ks results from scattering lIlSll . The function (fT2l) has been 
reformulated from its original version [@] to provide increased stability for scattering- 
free simulations involving extremely small opacities (i.e. fc^ = fc ^ 1), as is typical in 
cosmology applications. 

In the Poisson equation for the gravitational potential the baryonic gas is coupled 
to coUisionless dark matter p^/,,, and the cosmic mean density (p) solely through their 
self-consistent gravitational field. Here, g provides the gravitational constant, and the 
dark matter density is evolved using the Particle-Mesh method described in II 19ll2d . 12111. 



ALGORITHM 

Instead of working with the equations directly in CGS units, we first normalize the 
system to render the values tractable for floating point computation. To this end, we 
define the scaled units 

x = x/u^, g = g/ug, t = t/ut, (13) 

where the constants Ux, Ut and Ug correspond to the typical magnitudes of length, 
time and mass at the start of the simulation. We further define the density unit factor 
Ud — Ug/{ux)^ and velocity scaling factor m,, = Ux/ut- We note that due to our use of 
comoving length, these constants are all redshift-independent. The proper length values 
at any point in the simulation are therefore given by 

Voper =XUxa{t). (14) 

With these unit scalings, we define the normalized variables 

Ph = p/ud, Vfe = Vfe/Mv, e = e/ul, (15) 
E = E/ (udul) , Hi = rii/ud, 4> = ^/"v- 

The proper densities may be obtained from the comoving densities through the formulae 

Ue 



■Eproper — E /ci (f) — E 



a\t) 

Un 



n/,proper = n//a^(f) = (16) 

3/ ^ ~ 



Pi,proper = Pb/ « (0 = Pb 



a\t) 



With these rescaled variables, we rewrite our equations ([B-® as the normalized system 

dfpb + -yh-^Pb = --pb^-yb. (IV) 

a a 

dm + - (vfe ■ V) = -Ub - -^yp - -vf (18) 

1 2^ 1 1 

dfe + -Vh-Ve = — -e - —V ■ (p\b) - -Vb-V^ + G- A (19) 
a a apb a 

dfni + -V ■ (n,-v/,) = a/,;n,n; - n.-ff , / = 1, . . . , (20) 
a 

diE + -V- (E\b) = V - ( DVE)--E + 4717) -ckE, (21) 

a ^ ^ a 

V^ = ^{Pb + Pdm-{p)). (22) 

Here, d now refers to the derivative ^. For clarity of notation, all subsequent variables 
are shown without the ~ superscript, although all solver algorithms operate on the 
normalized variables. 



Operator-Split Hydrodynamics with Radiative Feedback 

We solve the coupled system (fT7l)-(l22l) using an operator- split framework, wherein we 
solve sub-components of the system one at a time, feeding the results of each sub-system 
into the remaining parts. In this approach, a time step is taken with the steps: 

(i) Project the dark matter particles onto the finite-volume mesh to generate the dark 
matter density field p^,,,. 

(ii) Solve for the gravitational potential (j) using equation (|22l) . 

(iii) Advect the dark matter particles with the Particle-Mesh method [ T9Ll20L l21]. 



(iv) Evolve the hydrodynamics equations (fT7])-(fT9l) with a high-order, explicit-time 
upwind method. In this step, the velocity field v/, advects both the chemical number 
densities n, and radiation energy density E. 

(v) Using a high-order implicit method, solve a coupled reaction-diffusion system (fT9l) - 
(|2TI) to obtain the updated number densities n,, radiation E and gas energy e. 



The equation (1191 ) is involved in both steps (iv) and (v) above. To do this, we split the gas 
energy into two parts, e = eii + ec, where e/, results from the hydrodynamic evolution of 
the system (step (iv)), and e^- is the gas energy correction that results from couplings with 
radiation and chemistry (step (v)). With this splitting, the hydrodynamic solver used in 



step (iv) of the algorithm solves the system of equations 



dtpb + -\b-^Pb = --pb^ -"Vb, (23) 
a a 

dtyb + - (Vfc ■ V) vz, = - -yb - — - - V0, (24) 
a a apb a 

1 2^ 1 1 

dteh + -yb ■ Ve/, = eh V • (jcv^) - -v^ ■ V0 (25) 

a a apb a 

<9,n, + ^V-(n,v^) =0, (26) 

<9,£ + ^V-(£vfc)=0, (27) 

using the Piecewise Parabolic Method (PPM) [f22|l . on a regular finite-volume 
spatial grid. This solve evolves {p'^,y'l,e'\x\!- ,E") to the time-updated variables 

(p^'+\vj'+\e)^+^) and the advected variables (n*,^*), and is implemented in the 



community astrophysics code Enzo [|21LI23L I24|] . 



Step (v) then solves the coupled system, 

2d 

d,ec = ec+G-A, (28) 

a 

d,ni = aijiieii j - HiVf , (29) 
djE = V ■{DVE)-m-E + AKT]-cKE, (30) 

using a fully implicit nonlinear solution approach to evolve the advected variables 
{e",rL*,E*) to the time-evolved quantities (e"+\n"^^, £■"+'). To this end, we define 
the vector of unknowns U = {ec,rii,E)^, and write the nonlinear residual function 
/(f/) = (/„/n,,/£)^, where 

MU) =ec-S,{U), (31) 
/n,(f/) =n,-5n,(t/), /=1,...,A^„ (32) 

fEiU) =E-E"-Ate (v-{DVE)-^E + 47tri-ckE^ (33) 

- Af(l - 0) ( V ■ {D"VE") - -E"+47tri" - ck"E' 
\ a 

Here the functions Se{U) and S^iiU) provide the analytical solutions to an 0{At^) 
accurate approximation of the spatially-local ODE system (l28l)-(|29l) for a given value 



of E n25\\ . The residual equation (|33l) defines a standard two-level method for time 
integration of the equation (l30l) . Therefore this overall nonlinear residual defines an up- 
to-second order implicit time discretization of the coupled system (|28l) -(l30l): where the 
time-evolved state U"'^^ is found through solution of the problem f{U) = 0. 

To solve this nonlinear problem, we use a globalized Inexact Newton 's Method ^ 



271] . that iteratively proceeds toward the solution U"^^ through approximately solving 



a sequence of linearized problems J{Uk)Sk = —f{Ui^), where J(Uk) is the Jacobian of 
the nonlinear function /, evaluated at the current Newton iterate U]^. h full description 
of our solution algorithm is provided in To summarize this process, we solve 
these linear Newton systems through a Schur complement formulation, that reduces 
the coupled linear system to a sequence of simpler sub-systems, culminating in an 
update to a modified radiation equation. This Schur complement subsystem is solved 
usinga multigrid-preconditioned conjugate gradient method from the HYPRE library 



We measure convergence of the Newton iteration with the RMS norm 



2 X 1/2 



where N{Ns + 2) is the number of unknowns in w, since this norm does not grow 
artificially larger with mesh refinement. 



VERIFICATION TESTS 

The model equations and solution algorithm in this paper have been rigorously tested 
on a myriad of problems, ranging from pure radiation transport, to interacting radiation 
and hydrodynamics, to dynamic radiation-hydrodynamics with chemical ionization 
In lieu of reiterating those tests here, we demonstrate the approach on a single problem 
before moving on to our target application. 

We consider a verification problem that performs isothermal ionization of a static (i.e., 
no fluid motions other than Hubble expansion) neutral hydrogen region, within a cosmo- 



logically expanding universe. The problem is originally due to Shapiro & Giroux 113111 . 
and exercises the radiation transfer, cosmology and chemical ionization components of 
the coupled solver. The physics of interest in this example is the expansion of an ion- 
ized hydrogen region in a uniform gas around a single monochromatic source, emitting 
A^y = 5 X 10^^ photons per second at the ionization frequency of hydrogen {hv = 13.6 
eV). Given the initially-neutral hydrogen region and strength of the ionizing source, the 
ionization region expands rapidly at first, with the I-front approaching the equilibrium 
position where ionizations and recombinations balance, referred to as the Stromgren ra- 
dius. However, due to cosmological expansion, this equilibrium radius begins to increase 
much faster than the I-front can propagate. The analytical formula for the location of the 
Stromgren radius as a function of time is 



rs{t) 



3Ny 



1/3 

(35) 



where the proper hydrogen number density lifj decreases due to cosmological expansion 
by a factor of a^^{t). Here Ub ~ 2.59 x 10^^^ cm-^/s is the case-B hydrogen recombi- 
nation coefficient. If we define A = agn//^o/^o/(l +zo), where n//.o is the hydrogen 
number density at the initial redshift zq, we may calculate the analytical I-front radius at 



TABLE 1. Cosmological parameters for the 
verification tests. See text for descriptions. 



qo ZQ Lq Hq n,„ D.A ilft 



0.5 4 80kpc 0.5 1.0 0.2 

0.05 4 60kpc 1.0 0.1 0.1 

0.5 10 36kpc 0.5 1.0 0.2 

0.05 10 27kpc 1.0 0.1 0.1 



any point in time as 



ri{t) = rs.o(^Xe-'^'^ J"' e^'^^l -2qo + 2qo{l +zo)/d]-'/^ da^ , (36) 

T{a) = X [6qlil+zo?r\F{a)-F{l)], (37) 
F (a) = [2 - 4qo - 2^o ( 1 + Zo) /a] [1 - 2qo + 2qo{l+ zo) /a] ^'^ , (38) 

and where qi^ is the cosmological deceleration parameter. 

We perform four of the tests provided in the original paper [Isil]: = {0.5,0.05}, 
and Zo = {4, 10}. These correspond to the cosmological parameters found in Table [B 
Here, Lq is the initial box size and Hq is the Hubble constant. The values Q.m, and 
are the contributions to the gas energy density at z = due to non-relativistic matter, 
the cosmological constant, and baryonic matter, respectively. These two deceleration 
parameters result in slightly different functions for the expansion coefficient a. For 
qQ = 0.05 we compute a{t) using equations (13-3) and (13-10) from [|32|]. while for 
qo = 0.5 we compute a{t) = (1 +z{t))^^. We begin all problems with an initial radiation 
energy density of E = 10^^^ erg cm^-^ and an initial ionization fraction nHii/n-H.o = 0. 
The initial density is chosen as p^ o = 1.175 x 10^^^ g cm^-^ for qo = 0.5, or p/, o = 
2.35 X 10^^^ g cm^^ for qo = 0.05. All simulations are run from the initial redshift 
Zo to z = 0. The ionization source is located in the lower comer of the box. We use 
reflecting boundary conditions at the lower boundaries and outflow conditions at the 
upper boundaries in each direction. The implicit solver used a convergence norm of 
p = 2, time step parameter of = 0.5 1 , a desired temporal solution accuracy Ttoi = 0.001 
and inexactness parameter 5^ = lO^^^\\f{Uk)\\ (see dH] for further explanation of these 
parameters). 

In Figure \T\ we plot the scaled, spherically-averaged I-front position with respect to 



where 



scaled redshift for each of the four tests (with axes identical to [|3 U\ . Figure la), as well 
as the zoomed-in version for the zo = 4 tests along with their analytical solutions; all 
of these tests used a 128^ spatial mesh. In Figure [21 we plot the error in the computed 
I-front radius as we varied the spatial mesh size for the two cases {qo.Zq) = (0.5,4) and 
(0.05,4). The accuracy in the computed radius improves with mesh refinement. 



l-front radius vs scaled redshift |_front radius vs scaled redshift, zO=4 




FIGURE 1. Left: I-front radius vs. scaled redshift for the four tests. Right: I-front radius vs scaled 
redshift for the zo~4- tests; analytical solution values are shown with open squares. 




FIGURE 2. Convergence of I-front radius vs. scaled redshift for the cases = {0.5,0.05} and zo = 4 
as the mesh is refined: spatial meshes shown are 16^ (blue dotted), 32^ (green dashed), 64^ (red solid), 
and 128^ (black dot-dashed). 



APPLICATION TO COSMOLOGICAL REIONIZATION 

To illustrate the operation of the combined cosmological radiation hydrodynamics plus 
ionization kinetics code, we simulate hydrogen reionization due to stellar sources in 
a small cosmological volume. This is a low resolution functionality test only to show 
that the two halves of the code are coupled correctly; scientific predictions will require 
considerably higher resolution and larger boxes [|33l1 . 

We simulate a ACDM cosmological model with the following parameters: = 0.7, 
Q.h — 0.04, Q.CDM = 0.26, h = 0.7, ag = 0.9, where these are, respectively, the fraction 
of the closure density in vacuum energy, baryons, and cold dark matter, the Hubble 
constant in units of 100 km/s/Mpc, and the variance of matter fluctuations in 8 Mpc^^ 



spheres, all measured at the present epoch (z=0). A Gaussian random field is initialized 
at z=99 using an initial power spectrum following [34]. The simulation was run in an 
8 Mpc comoving box using a 64^ uniform grid and dark matter particles with periodic 
boundary conditions. 

The emissivity rj was computed using a modified version of the star forma- 
tion/feedback recipe of [[l7|] . in which the conditions for star formation within a 
computational cell require that the local V/, have negative divergence (i.e baryons are 
contracting), the radiative cooling time is smaller than the dynamical time, and that ph 
is greater than some threshold (without checking for the Jeans mass). If these criteria 
are met a star particle is created which represents an ensemble of stars and becomes a 
source of emissivity for the radiation solver. Many particles are created over time and 
there may be multiple star particles in a cell. The emissivity of a cell rj is computed as 
follows: 



1 



^ = ^L^t^^/ 'nsF{t)c^dt (39) 

where the sum is over all the star particles in the cell, msF{t) is the star formation rate, 
which is an assumed analytic function of time for each particle, c is the speed of light, 
Af is the timestep, and euy an efficiency factor that depends on a number of hidden 
parameters including the initial mass function of the star cluster, the stellar spectral 
energy distribution, and the ionizing photon escape fraction. For simplicity, we used the 



upper value from [|16l] in its place. 



Snapshots from the simulation are created by the analysis tool yt [|35ll and shown 
in Table [2l Here projections through the three dimensional volume are shown for 
four redshifts z = 4.35, 2.55, 1.99 and and three physical quantities: baryonic den- 
sity (logio(pfc/Pav^)), ionization fraction (logio(n/////nH)) and temperature (log^jT) in 
Kelvin. 

The first star particle was created atz = 5.58. At this point, the initially homogeneous 
Intergalactic Medium (IGM) had formed filamentry structures as a result of dark matter 
clumps. By the first snapshot at z = 4.35, the effects of the star can clearly be seen in 
the higher ionization fraction in the lower right corner region around the star. This same 
region is evident by a brighter peak in the density and temperature. 

By z = 2.55, multiple sources have formed and are contributing to the ionizing 
radiation. The ionization fronts have also clearly propagated through significantly more 
of the domain. Although the ionization fronts are converging, there are still small pockets 
where the IGM remains neutral. By inspection, the majority of the IGM is at around 10'^ 
K, consistent with expectation. The bright peaks in temperature mark a region with 4 
hot stars in close proximity, but the area of neutral IGM remained cooler. 

A short while later (on the cosmological scale) at z = 1 .99, there has been little change 
in the density structure, but the ionization fronts have passed one other, overlapping 
the ionization region. At this point the universe is becoming transparent to ionizing 
radiation. Although in this simulation reionization has finished much later than observed 
[0], we note that the redshifts at which stars are made by this star maker recipe are 
heavily dependent on the box size and spatial resolution of the simulation. A bigger box 
size with the same spatial resolution will likely create stars at a much earlier redshift. 



TABLE 2. Cosmological ionization snapshots; 2D data results from averaging through one direction. 
The rows correspond to times z = 4.35, 2.55, 1.99, 0; the columns show baryonic density, ionization 
fraction and temperature. 
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Meanwhile, most of the computational volume has reached the same temperature, aside 
from the local hot spots around the stars. 

Finally at z = 0, the matter has nearly all coalesced due to the gravitational potential 
of the high density peaks. This results in large voids of underdensity regions in the IGM, 
and at the same time multiple bright spots are converging towards each other. By this 
time we see that the universe has been completely ionized (the data shows very small 
specks where some HI exists, which may be attributed to recombination). Furthermore, 
the temperature plot shows that most of the IGM is actually at a lower temperature than 
earlier at reionization. This is due to the adiabatic expansion of the universe, causing 
regions far away from the sources of radiation to cool. The brighter temperature region 



is also expanding. This is not due to the photo-ionization of the IGM as before, but is 
instead due to collisional heating from infall onto the massive dark matter halo, shock 
heating the regions around it. In a color plot, it can be seen that the shock front has a 
higher temperature than the relaxation area behind the front. 

Our box is far too small, and the spatial resolution too low, to describe the z = 
structure of the universe accurately. Indeed, density fluctuations with wavelengths 
comparable to the box size are going nonlinear at z = 0, making our solution highly 
inaccurate. The sole purpose of continuing the calculation was to test the long term 
stability of our implicit algorithm. It passed the test. 

CONCLUSIONS AND OUTLOOK 

The combined code appears to be working as expected and is stable for long execu- 
tions. Radiation from star formation fully ionizes the volume, however the redshift of 
reionization is delayed due to the low spatial resolution which underestimates the star 
formtion rate. Higher resolution runs and larger box sizes are planned in the near fu- 
ture. As discussed in [U] our radiation solver is optimally scalable with respect to the 
number of radiation sources, the number of grid points, and the number of processors. 
Moreover, the timesteps for the radiation-ionization kinetics portion of the calculation 
is independent of resolution because of the implicit time differencing. This is not the 
case for explicit cosmological dynamics, which means that at some grid size the radia- 
tion portion of the calculation will cease to dominate the runtime cost. We have not yet 
determined where this crossover occurs, but are investigating the matter. 

Several extensions of the method are under development. The first is multigroup 
FLD for a more accurate representation of the transport of hard UV and X-ray photons 
and helium ionization. A second is replacing the FLD ansatz with the variable tensor 
Eddington factor method used in [|36l] . This will improve the angular description of the 
radiation field and allow for shadowing effects. 

Finally, there is extending the radiation-ionization kinetics solver to adaptive mesh 
refinement. There are two components in this solver that depend on the spatial mesh. 
The first of these is the solver for the Schur complement subsystem. The part of the 
current solver for this component that currently depends on a uniform spatial mesh is the 
geometric multigrid solver that is used to precondition the conjugate gradient iteration. 
In extending the approach outlined here to spatially adaptive meshes, this geometric 
multigrid solver may be replaced with a Fast Adaptive Composite (FAC) method that 
understands the overall composite mesh that is formed out of a nested hierarchy of 
uniform grids of different spatial resolution. 

The second component that depends on a uniform spatial mesh is the rather straight- 
forward operator- splitting approach coupling the explicit and implicit sub-solvers. Due 
to the mesh-dependent CFL stability restriction, the explicit solvers employ time subcy- 
cling on the composite mesh, wherein more highly refined grids use smaller time steps 
than their larger counterparts, synchronizing with one another only at the time step of the 
coarsest grid. The implicit solver, however, naturally couples all of these levels together 
at once. Therefore in extending these solvers to AMR, we plan to examine the proper 
operator- splitting strategy for coupling these solvers together, attempting to balance a 



need for accuracy and consistency (use a full implicit solve every subcycled time step) 
with a need for efficiency (use a full implicit solve only at the coarsest grid time step). 
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